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Abstract 

We have performed Monte Carlo simulations to determine current-voltage 

characteristics of two different vortex glass models in two dimensions. The 

results confirm the conclusions of earlier studies that there is a transition at 

T = 0. In addition we find that, as T — > 0, the linear resistance vanishes 

exponentially, and the current scale, J n i, where non-linearities appear in the 

I—V characteristics varies roughly as T 3 , quite different from the predictions of 

conventional flux creep theory, J n i ~ T. The results for the two models agree 

quite well with each other, and also agree fairly well with recent experiments 

on very thin films of YBCO. 
PACS: 74.40.+k, 74.60.-w, 74.76.-w 
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I. INTRODUCTION 



Fluctuation effects are much larger in high T c superconductors than in conventional 
materials Q, so there has been considerable effort, recently, to go beyond the traditional 
mean field approximation in describing the superconducting phase transition. The effects of 
fluctuations, and also of disorder, are particularly strong in a magnetic field. In fact much 
of the H — T phase diagram of high T c materials is occupied by a "vortex liquid" regime in 
which the resistance has dropped because superconducting short range order has formed, so 
flux lines exist locally, but the resistance is not yet zero because the flux lines move under 
the action of a Lorentz force due to the current, and hence give rise to a voltage 0. An 
important question is whether, at lower T, the flux lines can be collectively pinned by defects 
sufficiently strongly that they have no linear response to the Lorentz force, which implies a 
vanishing linear resistance. In the presence of disorder, there is no long range order in the 
arrangement of vortices ||. Nonetheless, it is argued j|] that a transition to a state with 
vanishing linear resistance can occur. Such a state is called the vortex glass. 

While there is theoretical |5||| and experimental evidence for a finite vortex glass 
transition temperature, T c , in bulk superconductors, several simulations |]8]|| have clearly 
shown that T c = in two-dimensional systems, and rigorous analytic arguments |9[] have 
established that there is no vortex glass order at finite T in 2-d. 

Recent work [|10| has shown that inclusion of gauge field fluctuations (ie, screening) 
changes the universality class, but one still has T c = in two dimensions. Even though 
T c = 0, the correlation length diverges as T — > 0, which leads to observable consequences at 
finite temperatures, as we shall discuss in detail here. In contrast, a different simulation of 
the non-linear IV characteristic of the d=2 gauge glass model found evidence for a finite T c . 



11]. However behavior consistent with T c = has recently been observed [12] in experiments 



on very thin, 16A, films of YBCO. 



Although the experiments |T2[ on the 2-d samples are in quite good agreement with the- 
ory, they measure different quantities from what has been calculated so far, and a scaling 
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hypothesis is needed to make a comparison. The experiments determined current-voltage 
(I-V) characteristics while the simulations investigated the size-dependence of the rigidity 
of the system with respect to a twist in the phase of the condensate. It therefore seems 
worthwhile to make a direct comparison by calculating I-V characteristics from the simula- 
tions on 2-d systems. As an additional benefit, we study two models which are somewhat 
different microscopically: the gauge glass, defined in Eq. (fj) below, which has effectively 
random fluxes penetrating the sheet, and a more realistic model, defined in Eq. ([14]) below, 
which has a net uniform field penetrating the sheet and a random pinning potential for 
the vortices. We find that universal properties are the same for these two models. Since 
the gauge glass has been extensively used for numerical studies of the vortex glass transi- 
tion it is reassuring that it gives the same results as a model with a net field, at least in 
two-dimensions . 



II. THE MODELS 

The first model that we study is the gauge glass |5]|||8|], whose Hamiltonian is 

^ g g = - Z cos (^ _ fa ~ A ij) • (!) 

The phase of the condensate, fa, is defined on each site, i, of a square lattice, with N = L 2 
sites. The sum is over all nearest neighbor pairs on the lattice. The effects of the external 
magnetic field and disorder are represented by the quenched vector potentials, A^, which 
are taken to be independent random variables with a uniform distribution between and 
2tt. 

To compute the I-V characteristics we need to incorporate dynamics into the model. 
The standard way of doing this is to view the model as a set of coupled Josephson junc- 
tions Hl3|-[T5|]. Josephson's and Kirchoff's equations for the current are then 



A? = TT + !c sin (0i - <t>j - Aij) + Vij(t) (2) 

Kq 



^ "j lij liiext j (4) 
5 

where i and j are nearest neighbor pairs. Eq. (0) expresses the sum of the current from 
site i to neighboring site j as the sum of a resistive current given by Vij/R , a Josephson 
current, and a Langevin current noise source f]ij{t). I c is the maximum Josephson current 
of the nearest neighbor pair, and Rq is the shunt resistance of the pair. The thermal noise 
has a Gaussian distribution with the following properties: 

=0 , (5) 

(Vij(t)Vki{t r )) = 2k ^ T 5ij,ki5{t - t') , (6) 
K 

which ensures that the system comes to thermal equilibrium at temperature T. Eq. @ 
is the Josephson relation connecting the voltage with the time derivative of the phase 
difference, 0j — <pj. Eq. (J|) is Kirchoff's law expressing current conservation at site i. Ii :ex t is 
the external current at site i. This is zero except for sites on the top row where an external 
current J = I / L, is fed in, and sites on the bottom row where the same current is extracted. 
The total current through the sample is then I. The average voltage across the system, V, 
is given by 

v = ~ E E !<*-*> P) 

^ ze iebottomjetop ab 

where the brackets, (...), denote a time average. The average electric field is then obviously 
given by E = V/L. For this model we work in units where h/(2e) = Ro = I c = 1. 
Throughout the paper we also set Boltzmann's constant to be unity. The equations of 
motion are integrated using a first order approximation with a time step of St = 0.05r, 
where r = h/(2eRoI c ) is the basic unit of time (which is set equal to unity). The scale of 
r is the typical time for a neighboring pair of sites to accumulate a relative phase of order 
unity. 

We have also studied an equivalent form for the gauge glass written in terms of vortices. 
To obtain this we replace the cosine by the periodic Gaussian (Villain) function |TjJ, and 
perform standard manipulations |16|JT7|], obtaining 



K = ~\ E("* - b ^ - 3){nj ~ h) , (8) 

where the {rii} run over all integer values, subject to the "charge neutrality" constraint 
Yji n i — 0, and are interpreted as the vortex "charges", and G(i—j) is the vortex interaction, 

G(i - j) = (-) V 1 " 6Xp[ik " (ri " rj)] (9) 
1 J) \LJ ^ A-2cosk x -2cosk y ' y ' 

At large distance, G(i — j) — > 27r log | r i — rj | . The vortices sit on the sites of the dual lattice, 
which lie in the centers of the squares of the original lattice. The magnetic fluxes 6, are the 
lattice curl of the vector potential and are given by (l/2iv) times the directed sum of the 
quenched vector potentials on the links of the original lattice which surround the site i of 
the dual lattice. 

We study the I-V characteristics of this vortex model by using Monte Carlo dynamics. 
That is, we equate Monte Carlo 'time' and real time, an approximation which is expected 
to be good in the limit of overdamped dynamics and which has proven reasonable in other 
simulations. |18|] Choosing a nearest neighbor pair, i,j at random, we try to increase rii 
by 1 and decrease rij by 1, thus transferring a unit vortex from j to i. If the change in 
energy is AE, the move is accepted with the probability 1/(1 + exp((3AE)) appropriate 
to the "heatbath" algorithm. An applied current density J gives a Lorentz force of Jh/2e 
on a unit vortex. This can be incorporated into the Monte Carlo moves ]Tj| by adding to 
AE an amount Jh/(2e) if the vortex moves in the opposite direction to the Lorentz force, 
subtracting this amount if it moves in the same direction, and making no change in AE if 
it moves in a perpendicular direction. Biasing the moves in this way takes the system out 
of equilibrium and causes a net flux of vortices in a direction perpendicular to the current. 
This then generates a voltage V, where 

v = y e (i v (t)) , (10) 



with 



the vortex current. Here t denotes a Monte Carlo "time" (incremented by At after each 
attempted move), and AQY(t) = 1 if a vortex at site % moves one lattice spacing in the 
direction of the Lorentz force at time t, AQj(t) = —1 if the vortex moves in the direction 
opposite to the Lorentz force, and AQj (t) = otherwise. We set At = 1/4N so that an 
attempt is made to move each vortex once in each direction, on average, per unit time. We 
shall use units where h/2e = 1 when dealing with vortex models. 

The linear resistance can also be obtained from the Kubo formula for fluctuations in 
the voltage in the absence of any net current. This is exact for discrete time Monte Carlo 

i.e. 



dynamics provided the sum over time is made symmetrical about t — [19 



1 oo 

R Vm = _ £ At (V(f)V(0)) , (12) 

/J t=-oo 

which, in our units, can be expressed in terms of the vortex current as 

1 oo 

Rnn = iff, £ A * (^(*K V (0)> ■ (13) 

Ll t=-oo 

Using Monte Carlo dynamics should be a good approximation near a critical point, where 
the vortex motion is slow and overdamped. However, because of discretization of time, and 
the fact that the fastest that a vortex can move is one lattice spacing per time step, it is 
not very satisfactory for large currents or high temperatures, where the vortex motion is 
ballistic. For example, at high temperatures with no bias current, a vortex moves in the ± 
x direction with probability, 1/4, so, from Eq. (0), the resistance is given by Ru n = 1/(2T) 
and tends to zero, which is unphysical. 

The gauge glass represents a system with random fluxes penetrating the film but with 
zero net field. It is a convenient model to study but it should be verified that it is in the 
same universality class as experimental systems which have a net uniform field penetrating 
the film and a random pinning potential for the vortices. We have therefore also studied a 
random pinning model with the following Hamiltonian, 

^p = "!E ri « G ( l "jK'~E t, W n W 2 > ( 14 ) 



where the n{i) are restricted to the values 0, ±1, G(i — j) is given by Eq. @, and v(i) is 
a random pinning potential, uniformly distributed in the interval —A < v(i) < A. We set 
A = 7r and fix the net filling, / (= (1/N) Y^i n i )■> which effectively determines the magnetic 
field, to be 1/4. We obtain I-V characteristics from Monte Carlo dynamics. At each Monte 
Carlo time we try to insert a +1,-1 pair at a randomly chosen pair of sites % and j. The 
analysis is then precisely the same as described above for the vortex representation of the 
gauge glass. 



III. SCALING THEORY 

To analyze the results it is necessary to understand how the I-V characteristics vary in the 
vicinity of a second order phase transition. A detailed scaling theory has been developed |J 
and we now summarize the results of this for the case of a zero temperature transition, 
where the correlation length diverges as 

Z~T-» , (15) 

and the relaxation time, r, also diverges. Normally one defines a dynamic exponent, z, 
by t ~ £ 2 , but since the transition is at T = 0, the relaxation has an activated form and 
diverges exponentially as T — > 0. Formally this corresponds to z = oo. 

The vector potential, A, enters the Hamiltonian in Eq.(P in the dimensionless form 
Aij ~ J 3 'A(r) • dr. Therefore A scales as l/£. The electric field is given by E = —d t A 
and so scales as l/(£r). Now JE is the energy dissipated per unit volume per unit time. 
The natural unit of energy is ksT so E • J scales like kgT / (t^r) and hence J scales like 
knT j^ d ~ x . It is important to keep the factors of T because T c = 0. Combining these results 
we obtain (for ks — 1) 

^-*m • 

where g is a scaling function. In d = 2 with T c = this becomes 
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Tt i = ^) ■ < 17 > 

From Eq. (pj]) one sees that the characteristic current scale, J n i, at which non-linear behavior 
sets in, varies with T as 

Jnl~T 1+ » . (18) 

Now the linear resistance is just 

R hn = lim ^ , (19) 
and g(0) must be a constant, which we take to be unity, so we can write 

db = 9 (r^) • (20) 

and 

TR lin = - = A exp (-AE(T)fT) , (21) 
r 

where AE(T) is the typical barrier that a vortex has to cross to move a distance £. One 
conventionally defines a barrier exponent ip by AE ~ £^ ~ j 1 -^ j n terms of which 

Ti? lin ~ exp (-C/T 1+ ^) . (22) 

It has been suggested |I| that ip may be zero in d — 2, leading to barriers which are either 
finite as T — > or diverge logarithmically. In the latter case the linear resistance would vary 
as exp (— C\ log(T)| /i /T), where \x is another exponent. 

In Eq. the linear resistance is seen to vanish exponentially fast as T — > 0. In 

these circumstances it is generally more difficult to estimate the form of possible power law 
prefactors than the form of the leading exponential dependence. We have incorporated a 
factor of T on the LHS of Eq. ( pT|) because it emerged naturally from the scaling ansatz. 
This factor of T also looks reasonable when compared with the Kubo formula, Eq. ([13]), since 
it is TRn n which is given by the voltage fluctuations. However, it is possible that additional 
factors of T could be present in the scaling region. 

8 



In a finite system, the I-V characteristics will also depend on the size of the system 
when the bulk correlation length, £, becomes comparable with L. According to finite size 
scaling [^UJ only the ratio L/^is important and so, from Eq. (|T5|), one can generalize Eq. ( |20| ) 
to 



E f J 
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L 1/u T) . (23) 



JR Vm * VTi+- 

The non-linear behavior in the finite-size regime, described by Eq. (^), is rather complicated 
because it involves a function of two variables. We have therefore studied non-linear behavior 
either in the range where finite size corrections are negligible or, having already determined 
v, choose L and T such that the second argument, l}l vr I ', is constant. 



IV. RESULTS 



In Fig. [1] we show some of the experimental data of Dekker et al. [12] on 16A films of 
YBCO. For small current densities, the data is flat, indicating Ohm's law, with a linear 
resistivity which decreases rapidly with decreasing temperature. As J is increased the data 
starts to curve upwards, indicating non-linear response. The current scale, J n i, at which 
non-linear behavior sets in, is also seen to decrease with decreasing temperature. Analyzing 
all their data, Dekker et al. find the linear resistivity varies as 

R Vm ~ exp(-C/TP) (24) 

with p ~ 0.6. This variation, which is less rapid than an Arrhenius form, p — 1, is difficult 
to understand from the classical models which we study here. It may therefore indicate that 
quantum tunneling of vortices is important at the lowest temperatures. On the other hand, 
J n i is found to vary with T like T 3 , which, from Eq. (|1^) implies v ~ 2, in agreement with 
the earlier simulations [^,0 which studied a classical model. 

We next discuss our numerical results for the I-V characteristics of the gauge glass in 
the vortex representation. Data for T times the linear resistance against 1/T is shown in 



9 



Fig. Q on a log-linear plot. A simple Arrhenius form, i.e. a temperature independent barrier 
height, AE, would correspond to a straight line. In fact for the largest size, the data is very 
close to a straight line indicating that the barrier exponent, ip, in Eq. (^) is zero or close 
to it. These results are consistent the suggestion |IJ of a logarithmically increasing barrier, 
but this weak dependence will be difficult to see on data which are over a modest range of 
temperature in finite size systems, 

Data for the non-linear response at finite J is shown in Fig. |[ The data for the smallest 
current density, J, was actually obtained for J = from the Kubo formula, Eq. (|13D . As 
in the experimental results in Fig. [I], one sees Ohm's law at small J (where each dataset is 
horizontal) with a linear resistance which decreases rapidly with temperature, as shown in 
more detail in Fig. ^|, but deviations from Ohms law occur at a scale, J n i, where the data start 
to curve upwards. One sees that J n \ decreases with decreasing temperature, as expected. 
From Fig. ^|, finite size effects appear quite small for L = 16 at T = 0.5 and 0.8. Assuming 
v ~ 2, the data for L = 32 in Fig. |2| at T = 0.35 should also not be significantly affected by 
finite size effects. We have therefore analyzed the data in Fig. [3] according to the expected 
result for bulk behavior, Eq. (^). The scaling plot is shown in Fig. [|. The data scales 
reasonably well with v = 1.8 which is in quite good agreement with other estimates 
The values of i?ii n in this plot were obtained from the Kubo formula. 

Next we discuss the data obtained for the phase representation of the gauge glass. Results 
for the linear resistance are shown in Fig. |5|. The data is consistent with an Arrhenius form 
for the largest sizes, as was found for the vortex representation (Fig. 0). Results for non- 
linear current voltage characteristics for the gauge glass model in the phase representation 
are shown in Fig. |6|. Unlike the vortex representation, which has discrete time and so a 
maximum vortex velocity, the dynamics of the phase representation uses continuous time 
|2~T| and so can sensibly be applied for large values of J (and also high-T). In Fig. |6| 
one clearly sees a "flux-flow" regime for large J, where the Lorentz force is sufficient to 
overcome pinning, and the only hindrance to vortex motion comes from friction. This leads 
to a resistance which is roughly independent of J and also only rather weakly dependent 
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on temperature. For small J, however, the vortices are pinned by defects and move only by 
activation over barriers. The linear resistance in this "flux-creep' region, is therefore much 
smaller than that observed for larger J and is also strongly temperature dependent. 

We next describe our results for the random pinning potential model in Eq. The 
linear resistance is shown in Fig. [7|. As for the gauge glass, the data for the larger sizes seem 
to be tending towards an Arrhenius behavior. 

The non-linear behavior of the random pinning potential model is shown in Fig. || As for 
the gauge glass results in Fig. ^ one sees a "flux-flow" regime for large J, and a "flux-creep" 
region at small J where the linear resistance is small and strongly temperature dependent. 
Deviations from the Ohm's law behavior in the flux creep regime occur at a current scale, 
J n i, which decreases as T decreases. For sufficiently large J, E will saturate, because there 
is a maximum vortex velocity when each vortex hops every step, and so the ratio E/J will 
decrease. This (unphysical) behavior, which is just visible in the figure at the largest values 
of J, is caused by discretization of time in the Monte Carlo simulations. This discretization 
is, however, not expected to affect universal critical properties near the T = critical point. 

A scaling plot of the non-linear I-V characteristics of the random pinning potential model 
is shown in Fig. Since there are finite-size effects within the range of accessible sizes, we 
assume that v ~ 2, and choose sizes and temperatures such that LT 2 is constant. In this 
way, the second argument in the scaling function in Eq. (|23|) is constant and E/(JRn n ) 
should only be a function of J/T 3 . The data is seen to scale very well over a wide range. 
These results provide strong evidence that the gauge glass and random pinning potential 
models are in the same universality class with a correlation length exponent at the T = 
transition of v ~ 2. 



Finally, in Fig. [L0[ we compare the non-linear current-voltage characteristics from the 
simulations on the random pinning potential and gauge glass models with the experimental 
results shown in Fig. [[]. The same value of v = 2 was used for all the data, and, for each 
system, a temperature scale, T , was adjusted to get the best scaling. The two sets of 
simulation data agree quite well with each other. While the gauge glass data were obtained 
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in the region where ^ < L, as deduced from the results for Ru n , the random pinning 
potential model data was in the finite size region and so the data were taken at fixed LT 2 ' 
(or equivalently at fixed, but not very small, £/L). The good agreement between the two sets 
of data indicates that finite size corrections are not very important for the I-V characteristics. 
There is, however, a discrepancy between the numerical results and experiment, for which 
we do not have an explanation. 

V. CONCLUSIONS 

We have studied the I-V characteristics of two models for vortex glass behavior in two 
dimensions. For the first model, the gauge glass, our results confirm earlier studies ||[8| which 
found a zero temperature transition with a correlation length exponent v ~ 2. This behavior 
has also been seen experimentally [jT^] on very thin films of YBCO, though the detailed form 
of the I-V scaling function is somewhat different between theory and experiment. We also 
find that the linear resistance varies with an Arrhenius form as T — > 0, indicating that the 
barrier exponent, ip, is either zero or close to zero. Experimentally the resistance appears 
to vanish at low T less rapidly than an Arrhenius form, which may indicate that quantum 
fluctuations of the vortices play a role in the YBCO films. The other model which we study, 
the random pinning potential model, is rather more realistic in that it describes a system 
with a uniform applied magnetic field perpendicular to the film, as opposed to the gauge 
glass which has random fluxes. Nonetheless, the random pinning potential model is found to 
be in the same universality class as the gauge glass, since they both have a zero temperature 
transition with v ~ 2 and an Arrhenius behavior for the linear resistance. Hence these 
two models are equivalent in d = 2. It can also be shown that in 6 — e dimensions, the 
gauge glass transition is in the same universality class as the vortex glass transition in a 
disordered model with a net magnetic field. It is therefore quite plausible that the two 
models have the same critical behavior in three dimensions, though we are not aware of a 
direct demonstration of this. 
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FIGURES 

FIG. 1. A plot of some of the experimental data, at a field of 0.5 T, on very thin, 16A, films 
of YBCO from Ref. The numbers denote the temperatures, in Kelvin, at which the different 
data sets were taken. 

FIG. 2. A plot of T times the linear resistance, on a logarithmic scale, against 1/T for different 
sizes for the gauge glass in the vortex representation. The data for the largest size, L = 32, is 
well approximated by a straight line, indicating a temperature independent barrier height, i.e. 
Arrhenius behavior. 

in 

FIG. 3. A logdog plot of the nondinear current voltage characteristics of the gauge glass model 
in the vortex representation. The data shown for the smallest current density, J, is actually for 
J = and was obtained from the Kubo formula, Eq. (|l3|). One sees Ohm's law at small J (where 
each dataset is horizontal) with a linear resistance which decreases rapidly with temperature, as 
shown in more detail in Fig. ^. Deviations from Ohm's law, where the data start to curve, occur 
at a current scale, J n i, which decreases with decreasing temperature. 

FIG. 4. A scaling plot of the data in Fig. assuming that finite-size corrections are small 
and hence that the scaling form expected for bulk behavior, Eq. (ptj|), is appropriate. The value 
v = 1.8 obtained from this fit, is is reasonable agreement with other estimates. 

FIG. 5. A plot of the linear resistance, on a logarithmic scale, against 1/T for different sizes 
for the gauge glass in the phase representation. The data for the largest size, L = 16, is fairly close 
to a straight line, indicating a temperature independent barrier height, i.e. Arrhenius behavior. 
This is similar to what was observed in the vortex representation, Fig. ||. 
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FIG. 6. Results for non-linear current voltage characteristics for the gauge glass model in the 
phase representation. One sees a "flux-flow" regime for large J, where the resistance is largely 
independent of J and T, and a "flux-creep" region at small J where the linear resistance is small 
and strongly temperature dependent. Deviations from the Ohm's law behavior seen at small J 
occur at a current scale, J n i, which decreases as T decreases. 

FIG. 7. A plot of the linear resistance, on a logarithmic scale, against 1/T for different sizes for 
the random pinning potential model in Eq. (|l4|). The data for the largest sizes seems to be tending 
to a straight line, indicating a temperature independent barrier height, i.e. Arrhenius behavior. 
This is similar to what was observed for the gauge glass, see Figs. | and |. 

FIG. 8. Results for non- linear current voltage characteristics for the random pinning potential 
model in Eq. (|14|). For each size, L, the temperature has been chosen so that LT 2 = 2, in order to 
keep the second argument of the finite-size scaling function in Eq. (|23|) constant. The net filling, 
/ (= (l/N)J2i n i )j is equal to 1/4. As for the gauge glass results presented in Fig. || one sees a 
"flux-flow" regime for large J, and a "flux-creep" region at small J where the linear resistance is 
small and strongly temperature dependent. Deviations from the Ohm's law behavior seen in the 
flux creep regime occur at a current scale, J n i, which decreases as T decreases. 

FIG. 9. A scaling plot of the non-linear current-voltage characteristics of the random pinning 
potential model. Since there are finite-size effects in this data, we assume that v ~ 2, and choose 
sizes and temperatures such that LT 2 is constant. In this way, the second argument in the scaling 
function in Eq. (23) is constant and E/{JRn n ) should only be a function of J/T 3 . The data is 



seen to scale very well to this expected form over a wide range. 
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FIG. 10. A scaling plot combining the non-linear current-voltage characteristics of both the 
experimental data in Fig. [l] (lines) and simulations (points). The data points with LT 2 = 2 were 
obtained for the random pinning potential model and the other two sets of data points, at T = 0.35 
and 0.50, were for obtained for the gauge glass model in the vortex representation. The value v = 2 
was used for all the data. The temperature scale was set by To = 1 (random pinning potential 
model), To = 1.15 (gauge glass), and To = 12K (experiment). The two sets of data from the 
simulations agree quite well with each other, but the experimental results lie lower than theory at 
large values of J/T 1+u . 
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